Spike propagation for spatially correlated inputs 
through noisy multilayer networks * 



Hideo Hasegawa ^ 

Department of Physics, Tokyo Gakugei University, Koganei, Tokyo 184-8501, Japan 

(February 2, 2008) 



Abstract 



Spike propagation for spatially correlated inputs in layered neural net- 
works has been investigated with the use of a semi-analytical dynamical 
mean-field approximation (DMA) theory recently proposed by the author [H. 
Hasegawa, Phys. Rev. E 67, 041903 (2003)]. Each layer of the network is 
assumed to consist of FitzHugh-Nagumo neurons which are coupled by feed- 
forward couplings. Applying single spikes to the network with input-time 
jitters whose root-mean-square (RMS) value and the spatial correlation are 
(7/ and sj, respectively, we have calculated the RMS value {aom) and the cor- 
relation {som) of jitters in output-firing times on each layer m. For all-to-all 
feedforward couplings, som gradually grows to a fairly large value as spikes 
propagate through the layer, even for inputs without the correlation. This 
shows that for the correlation to be in the range of observed value of 0.1-0.3, 
we have to take into account noises and more realistic feedforward couplings. 
Model calculations including local feedforward connections besides all-to-all 
feedforward couplings in multilayers subject to white noises, have shown that 
in a long multilayer, aom and som converge to fixed-point values which are 
determined by model parameters characterizing the multilayer architecture. 
Results of DMA calculations are in fairly good agreement with those of direct 
simulations although the computational time of the former is much smaller 
than that of the latter. 
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I. INTRODUCTION 



In living brains, information is carried by spikes which propagate from one cortical area 
to another area. It has been controversial how information is coded in spikes (for review, see 
[1]- [7]). One possibility is that information is coded in the number of spikes within a short 
time window {rate code) [8]. Indeed, firing activities of motor and sensory neurons have 
been reported to vary in response to applied stimuli. In an alternative temporal code, on the 
contrary, information is assumed to be carried by precise firing times of neurons [9]- [11]. One 
of problems in the rate coding is that a fairly long time of tens of milliseconds are required 
to read out the rate for typical firings rate of 10-100 Hz. However, human visual systems, for 
example, have been reported to quickly classify patterns within 250 ms despite the fact that 
at lest ten synaptic stages are involved from retina to the temporal brain: transmission times 
between the two successive stages are no more than 10 ms on the average [12]. A possible 
mechanism to speed up reading of firing rate may be to collect spikes of many independent 
neurons in a population {population code) [13] [14], where many parallel neurons perform 
the same task with the inefficient, high redundancy. On the other hand, one of problems 
in the temporal coding is that spikes are vulnerable to noise while the rate coding performs 
robustly but inefficiently. These issues on coding have been theoretically studied in a single 
neuron ensemble. It is not clear whether these conclusions may be applied to multilayer 
architectures relevant to cortical processing. 

Studies with the use of multiunit recordings of frontal cortex of monkeys have shown 
that a spatiotemporal pattern of highly synchronous firings can propagate through several 
tens of synaptic connections [15]. A simple model accounting for this phenomenon is a 
feedforward synfire chain first proposed by Abeles [15]. Since the synfire chain model was 
proposed, many studies have been made on its properties [16]- [24]. In the rate coding, 
neurons in each layer are expected to fire in uncorrelated manner with other neurons in the 
same layer. Neurons in a given layer are assumed to compute the average firing rate of the 
neurons in the previous layer in order to generate the output rate which is related to the 
input rate. In a feedforward network, however, this firing may propagate to the next layer 
in synchronous way [15], which is detrimental for the rate code. Diesmann, Gewaltig and 
Aertsen [18] have shown by simulations of integrate-and-fire (IF) neuron model that a pulse 
packet can propagate through the synfire chain if a packet satisfies the condition which is 
specified by the two parameters: one is the number of spikes in a pulse packet and the 
other is the root-niean-square (RMS) value of firing times in a pulse packet. The result of 
Diesmann et al. [18] has been confirmed by the method of Fokker-Planck equation [19]. 

It has been not clear whether feedforward networks support the rate-code or temporal- 
code hypothesis. Shadlen and Newsome [25,26] have claimed the feasibility of the rate 
code, adopting a model in which excitatory and inhibitory synaptic inputs are assumed to 
be balanced. Because of this balanced input, postsynaptic potentials fiuctuate around the 
resting potential, which yields random firings in output neurons. It has been shown that if 
each pair of output neurons shares less than 40 percent of input neurons, only a small degree 
of synchrony will be developed, which assures an feasibility of the rate code. The efficiency 
of the rate code transmission in unbalanced feedforward networks has been also studied [22]. 
Quite recently, however, it has been pointed out that in long feedforward networks, input 
firing rate cannot be transmitted reliably because the mean firings rate in a deep layer is 
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independent of input firing rate [23]. 

Studies on feedforward networks have been so far made mostly by direct simulations for 
networks described by the simplest IF model. It is worthwhile to make a more detailed study 
on feedforward multilayers by employing more realistic neuron model with an analytical 
method besides simulations. In a previous paper [24] (referred to I hereafter), we have 
developed the semi-analytical dynamical mean-field approximation (DMA) theory as an 
efficient tool dealing with large-scale FitzHugh-Nagumo (FN) neuron ensembles subject to 
noises [27] [28], by extending the moment method [29]. Original 2A^-dimensional stochastic 
differential equations (DEs) for a A'^-unit FN neuron ensemble are transformed to iV(2A^+3)- 
dimensional deterministic DEs for means, variances and covariances of local and global 
variables. Recently DMA has been successfully applied to neuron ensembles described by 
the realistic Hodgkin-Huxley (HH) model [30] [31]. The FN neuron model adopted in I for 
a feedforward network is obtainable by a simplification of the HH model [27] [28] , and it is 
expected to be more realistic than IF model. We have investigated in I, the spike propagation 
through the network, taking no account of the spatial correlation. Experimentally, correlated 
firings have been observed in a variety of neurons [32]- [40]. It has been reported that the 
correlation coefficient between cells is about 0.12 in V5 of a rhesus monkey [32], 0.1-0.3 in 
human motor units of muscles [36], and about 0.3 in cat's lateral geniculate nucleus (LGN) 
[34] and in retinal ganglion cells of rabbits [39] . Theoretical studies on the input correlation 
have shown that it may yield a significant effect on the firing rate and the variability of 
outputs [41]- [48]. Calculations with the use of IF model have shown that the firing rate 
of outputs is increased with increasing the input correlation for low firing rates of inputs 
but is decreased for their high firing rates [41] [46] [48]. The variability of output spikes of 
IF model is an increasing function of the input correlation, whereas that of HH model is a 
decreasing function of the input correlation [47] . 

These studies have been made for a single neuron ensemble [41]- [48]. We expect that 
the spatial correlation plays an important role also in multilayer networks. Although some 
theoretical studies have investigated the cross-correlation of spike rates averaged over long 
times [23], there have been no calculations of the firing-time correlation in multilayers, as 
far as the author is concerned. We have developed, in I, a new method calculating the 
instantaneous synchronization ratio in neuron ensembles which is expressed in terms of 
variances of local and global variables [Eq. (58)]. As will be shown shortly, the calculated 
correlation in multilayers with all-to-all feedforward couplings subject to weak noises is 
developed to a fairly large value as spikes propagate, even for inputs without the correlation. 
In order that the correlation remains in the range of the observed value of 0.1-0.3 mentioned 
above [32]- [40], we have to take into account at least two factors: one is the more detailed 
connectivity in feedforward couplings besides the all-to-all coupling and the other is noises. 
As for the first issue, we have assumed, in this study, that our multilayer network includes, 
besides all-to-all couplings, local couplings in which each neuron in a given layer receives an 
input from one neuron in the preceding layer. All-to-all and local couplings are superimposed 
with fractions of p and 1 — p, respectively, where p denotes a parameter expressing a degree 
of all-to-all component in the total feedforward couplings. As for the second issue, several 
conceivable sources of noises have been reported: (i) cells in sensory neurons are exposed to 
noisy outer world, (ii) ion channels of the membrane of neurons and synaptic transmission 
by a release of synaptic vesicles are essentially stochastic, and (iii) synaptic inputs include 
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leaked currents from neighboring neurons. In this study, we have taken account of white 
noises which are independently added to all neurons. Applying spike inputs to the first layer 
of the network with the spatial correlation in input-time jitters, we have investigated the 
effect of the spatial correlation in multilayer with all-to-all and local feedforward couplings 
subject to independent noises. 

The paper is organized as follows: In Sec. II, we will discuss an adopted multilayer with 
feedforward couplings. By using DMA, the RMS value {(Tom) and the correlation (som) of 
jitters in output firing times at layer m are expressed as functions of the RMS value (cr/) 
and the correlation (sj) of jitters in input times. In Sec. Ill, some model calculations are 
reported of the correlated spike propagation through the multilayer by using DMA theory 
and direct simulations. The final Sec. IV is devoted to discussions and conclusions. 



II. LAYERED NETWORKS CONSISTING OF FN NEURONS 

A. Adopted model 

We have adopted M-layer neural networks in which each layer includes A^-unit FN neu- 
rons. Dynamics of a single FN neuron j (=1 to A^) in a given layer m (—1 to M) is described 
by nonlinear differential equations (DEs) given by 



Jt 
dymjjt) 
dt 



F[x^j{t)] -cy^,{t) + ltf(t) + + li]{t) + (1) 

bxmj{t) - dy,rij{t) + e, (2) 



with 



^!S\t)=(l^) j:G{Xmk{t)), (3) 

ltf{t) = (1 - 5m,) w, [(^) Y.G{xm-ik{t)) + {I - p)G{xm-i,{t))i (4) 
lf{t)^5miuH{t). (5) 

In Eqs. (l)-(5), F[x{t)] = 0.5 x{t) [x{t) - 0.1] [1 - x{t)], b = 0.015, c = 1.0, d = 0.003 and 
e = [24] [29], and Xmj and ymj denote the fast (voltage) and slow (recovery) variables, 
respectively, of a given neuron j in the layer m; I^^fi't) in Eq. (3) denotes the intra- layer 
couplings with the strength wi, the sigmoid function G{x) given by G{x) — l/[l-|-exp[— (x — 
0)/x]], the threshold 6 and the width x [49]; the first and second terms of Imfit) iii Eq. 
(4) stand for all-to-all and local couplings, respectively, with the inter-layer feed-forward 
couplings W2, p denoting the degree of common inputs to neuron j in the layer m from 
neurons in the preceding layer m — 1 [50]; jj'^'' in Eq. (5) denotes inputs applied to the 
first layer with magnitude of u and an arbitrary function of H{t) whose explicit form will 
be specified below [Eq. (8)]; the last term of Eq. (1), ^mj{t), expresses the Gaussian white 
noise given by 
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< Cmjit) > = 0, (6) 
< Cmjit) U{t') > = Sjk Smn Sit - t'), (7) 

where /3 denotes the magnitudes of noises and the bracket < • > expresses the expectation 
value. 

We will study effects of the spatial correlation in single spikes on the propagation of spike 
inputs. We adopt an external input Ij^'' in Eq. (5) with the alpha function, a{t): 

H{t) = a{t - tij) = [{t - tij)/rs] exp[l - {t - tij)/rs] e{t - hj), (8) 

where Ts stands for the synaptic time constant and Q{t) the Heaviside function given by 
G(t) = 1 for t > and otherwise. We assume that jitters in input times tjj in Eq. (8) 
obey the Gaussian distribution with means and variance given by 

< tij > =ti, (9) 
<StijStik> ^aj[Sjk+{l-Sjk)si], (10) 

where 5tij — tij —ti, ui and si denote RMS value and the spatial correlation, respectively, 
of input-time jitters. 

When an input spike given by Eqs. (5) and (8) is applied to the first layer, it may 
propagate through the multilayer in the propagating regime [18]. The firing time of a given 
neuron j in the layer m is defined as the time when the fast variable x^jit) crosses the 
threshold 9 from below: 

tomj = {t I x„,j{t) = 9; x^j > 0}. (11) 

Means, RMS value and the spatial correlation of jitters in output firing times on the layer 
m are given by 

tom = < tomj >, (12) 

(^Om ^ < Stomj >> (13) 



where Stomj — tomj — tom- We will calculate aom and som as functions of aj and sj for a 
set of model parameters by direct simulations and DMA theory, details of the latter being 
discussed in the following subsection. 

B. DMA theory 

1. Equations of motions 
As in I [24] , we first define the global variables for the layer m by 

X'^it) = ^ E (15) 
j 

^"^W = ^ E yn^M (16) 
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and their averages by 

/.r(t) = < x^{t) >, (17) 

/x^(t) = < Y"^it) > . (18) 

Next we define variances and covariances between local variables in the layers n and m, 
given by 



7iT(^) = 


1 

TV 


j 


(19) 


72T(^) = 


1 

TV 


j 


(20) 


7iT(^) = 


1 

TV 


j 


(21) 


72T(^) = 


1 

TV 


< ^Vnjit) SXmj{t) >, 


(22) 



and those between global variables in layers n and m, given by 

p^f (t) = < 5X"(t) (5X'"(t) >, (23) 

P2:T{t) = < SY^it) SY^'it) >, (24) 

(^) = < 3X^{t) SY^'it) >, (25) 

p^'^(t) = < 5y"(t) SX'^it) >, (26) 

where 5x„,(t) = a;^,(t) - /i™(t), <5|/„,(t) = y„^.(t) - /,™(t), SX^{t) =X'%t) -^^™(t) and 
5F™(t) = Y"^{t) — iJi^{t). It is noted that for n = m, we get 7^2"* = 7^1™ and p™2™ = ^2*1™- 
In deriving equations of motions, we have assumed small P and o"/, and the Gaussian 
distribution of state variables, as in I. The interlayer correlation between layers, which was 
neglected in I, has been taken into account within the nearest-layer approximation (NLA) 
in which the correlation beyond adjacent layers is neglected, as given by 

pZ''"" = 0. for ^ > 1 {k, a = 1, 2) (27) 

After some manipulations, we get the following DEs for m=l to M (arguement t is ne- 
glected) : 



IT 

dt 

7 m.m 



dt 

n,'t 

;,2 



7 m,m 
072,2 



dt 

n,'i 
,2 



dt 



+ fTltr - cf^2 + ^i9o + 5ml uh, + {l- 5m,) W2 g^-\ (28) 

hpJ^ - dn^ + e, (29) 

2(a™7D" - C7rr) + ^^idTCT + + 2^!?i, (30) 

2{b^^r-dj^r), (31) 

i.-,m,m , / m i\ m,in m,m . m/-m,m . -^m fQn\ 

07i,i + (« -c()7i_2 -C72,2 +Wi5(iCi,2 +Ai^2> (32) 
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= 2{a-pT:r - cpT,n + ^w^g^pTr + + 2^m' (33) 



m,m 



dp2,2 

dt 
dpi 



m,m 



2 



dt 



2{bpTf-dp-n, (34) 
= bpT,r + {a- - d)pT,r - cp7,r + wigTpTf + Y,-„ (35) 



with 



/ n.rn ^n,m / ]\t\ 

W - (1-1/iV) ' ^^^^ 

X™ = 5^1 u h Piit) + (1 - 5™i) ^2 ^7^"' [pPm + (1 - P) 7i:r''"], (37) 

X- = « ^2(t) + (1 - t^^2 ^7r"' b P™2"''"^ + (1 - P) Tm^'H, (38) 

= (5^1 /ii i?i(i) + (1 - (^^1) «;2 ' (39) 

= « ^1 Mt) + (1 - <^mi) ^r' (40) 

where a™ = /f + 3/3"7[';i, /^^ = (1/£!)FW(/x™), £ff = (l/£!)GW(/xf ) and he = 
{l/i\)d^ H{t)/dt^. In Eqs. '(37)-(40), P« and i?^ {k = 1,2) express contributions to the 
first layer [m — 1), obeying the following DEs (see the Appendix A): 

^ = a'P, - cP, + [NR, - P,] + a} u h,, (41) 

^ = 6Pi - dP2, (42) 

^ = a^R, - cR2 + w, gl i?i + + (1 - -^) s/] u (43) 
di?2 



dt 



bRi - dR2. (44) 



In Eqs. (37)- (40), ^'"^ and p^,^^'"^ express the interlayer correlation between layers m — 1 
and m, satisfying DEs given by 

7 m—l,m 

Jij-- = (a-1 + a™)7r:r''"^ - c(7r2-''"^ + i2,i'n + ^i(^r' + ^Dcr ''"^ 

+ ^2^rVPM ''"^-^ + (1 -P) 717''""']' (45) 

T2,2 7/ m— l.m , m— l,m\ oj m— l.m //icN 

= ^'(7i,2 + 72,1 ) - 2^72,2 , (46) 

m—l,m 

= &7r:r''"^ + (a""-' - dhTi''"^ - cj^i''"^ + w^gT-\^2-''"', (47) 

m—l,m 

2^ = 6717^'™ + (a™ - d)72":r'''" - C72":2-''™ + w^gTC^^''"' 

+ W2gr'\p pZ'^"^-' + (1 -p) 72T'''""']> (48) 
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+ w2gr'pT,i'''^-\ (49) 

m—l,m 

P2,2 7/ m— l,m , m— l,m\ r, r m— l,m /rr\\ 

= 2 + P2,i ) - 2rfp2,2 , (50) 

m—l,m 

Pi, 2 7 m— l,m I / m— 1 j\ m— l,m m—l,m , m— 1 m— l,m /criN 

- +(« -C?)Pl,2 -CP2,2 Pl,2 , (51) 

m—l,m 

P2,l 1 m—l,m I /„m i\ m—l.m ^m—l.m , „m„m— l.m , „m—l„m—l,m—l / rn\ 

= ' + (a - rf)p2,l - CP2,2 + WlQi P2,l + W29i P2,l ■ (52) 

Original 2A^M-dimensional deterministic DEs given by Eqs. (l)-(5) are transformed to N^q- 
dimensional deterministic DEs given by Eqs. (28)-(52) where Ngq = 12 + 16(M— 1). If there 
are no jitters (P^ — R^, — 0) or if the interlayer correlation is neglected (p"'^ = for m 7^ n), 
DMA leads to 8M-dimensional DEs as in I. We note that the contribution of input jitters 
to the first layer is proportional to aj in Eq. (41) whereas that to [1/N + (1 — 1/N)si]aj in 
Eq. (43). 

2. Quantities relevant to output firings 

Neuron Activity and Firing-Time Distribution 

We will show, in this subsection, that from iJi^{t), ^^^C{t) and p^f^{t) which are obtained 
from Eqs. (28)-(52), we may calculate the three important quantities relevant to output 
firings on the layer m: the activity of neurons (aom), the RMS value of output jitters (crom) 
and their spatial correlation (som)- 

The averaged distribution of the voltage variable Xmj{t) is described by the Gaussian 
distribution with the mean of fji^it) and the variance of ^^{^{t) [24] [29]. The probabihty 
Wom{t) when Xjnj{t) at t is above the threshold 9 is given by 

Wom{t)-l-^l^\'-^^\, (53) 



where il^{y) is the error function given by an integration from —00 to y of the normal 
distribution function (f){x): 

m = ^exp f-^^l . (54) 



/27r 

The neuron activity in the layer ni is given by 

aOm = Wom{t*Om)^ (55) 

which is unity when all neurons in the layer fire t — defined by /i7*(^om) — ^- The 
fraction of firings of neurons in the layer m is given by [24] 

Zomit) = ^ ^ 0(^-^) -|(-^) eiHn, (56) 
dt aom dt J^]^f^ 

with the RMS value of jitters of output spikes given by 
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7iT 



where /x™, /ii™" and 7™i™ are evaluated at t^^. Our cxom corresponds to cxo, RMS of firings 
times, of Diesmann, Gewaltig and Aertsen [18]. 

Synchronization ratio and Correlation of output firings 

The synchronization ratio Smit) in a given layer m is given by [24] 

Smit) (1 _ 1/^) (58) 

1 ^ ^^mj ^-^mk (59) 



N{N-1) . km ^f<5^l~><S^^' 



which is and 1 for completely asynchronous and synchronous states, respectively. Then 
the spatially-averaged correlation of output firing times in layer m defined by Eq. (14), is 
given by 

1 < StOmj StOmk > C ^-^* \ rRn\ 



where the relation given by Eq. (57) is adopted. 

Thus aom, CTom and som given by Eqs. (55), (57) and (60), respectively, are expressed 
in terms of /i^, 7™i™ and and they depend on model parameters of aj, sj, p, (3, wi, 

W2 and N. 



III. MODEL CALCULATIONS 
A. Effects of si 

In this study, we pay our attention to the response of multilayer networks to a single 
spike input of l'^'^\t) with tj = 100 in Eqs. (5) and (8). We have adopted the parameters 
of u — 0.10, 9 — 0.5, X — 0.1 and Tg — 5. Parameter values of ai, sj, p, /3, wi, W2, N 
and M will be explained shortly. The value of « = 0.10 has been chosen for a study of 
the response to a supra-threshold input, because the critical magnitude of u is Uc = 0.0435 
below which firings of neuron defined by Eq. (11) cannot take place for o"/ = = 0. Direct 
simulations have been performed by solving 2MN DEs given by Eqs. (l)-(5) with the use of 
the fourth-order Runge-Kutta method with a time step of 0.01 for hundred trials otherwise 
noticed. Correlated input times of tij given by Eqs. (9) and (10) have been generated by 
the Gaussian-distribution programs. DEs of DMA given by Eqs. (28)- (52) have been solved 
by using also the fourth-order Runge-Kutta method with a time step of 0.01. All calculated 
quantities are dimensionless. 

Raster in Fig. 1(a) shows firings of neurons of the first ten layers in a multilayer of N — 10 
and M = 20 for a typical set of parameters of aj = 1, sj = 0, p = 1, j3 — 0.01, ^1 = 
and W2 — 0.1, calculated by a direct simulation (a single trial). The ordinate expresses the 
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neuron index k defined hj k — 10(m— where m — 1 — 10 and j = 1 — 10. The uppermost 
chistcr denotes firings of ten neurons in the layer m = 1 and the bottom cluster those in 
the layer m = 10. When spikes are applied at t = 100, neurons are already randomized 
because noises have been apphed since t — 0. Firings occur with a delay of about 5 at each 
stage, and it takes about 48 for spikes to propagate from m = 1 to m = 10. When the noise 
intensity is increased to j3 = 0.02, fluctuations of firings due to inputs and spurious firings 
are increased, as shown in Fig. 1(b). Figures 2(a) and 2(b) show time courses of /i7^(t) and 
Sjn{t) of the first ten layers for the same set of parameters as in Fig. 1(a): solid curves 
denote the results of DMA theory and dashed curves those of direct simulations. Figure 
2(a) shows that a spike propagates from m = 1 to m = 10. Results of n^{t) of DMA are in 
good agreement with those of direct simulations; the former is not distinguishable from the 
latter. The synchronization ratio of Sm{t) shown in Fig. 2(b) is zero at m = 1 because of 
the vanishing input correlation s/ = 0. Nevertheless, Sm{t) after receiving inputs gradually 
become large as a spike propagates through the layer. This development in the synchrony 
is more clearly realized in Fig. 3(a), where large open and filled circles, respectively, show 
the m dependence of the correlation of Som = Sm{tom) S/ = and p = 1 calculated by 
direct simulation (dashed curve) and DMA (solid curve). We note that although som = at 
m = 1 for s/ = 0, it is rapidly increased and saturates with a value of about 0.71 (0.61) in 
direct simulations (DMA calculation) at m = 20. On the contrary, open and filled squares, 
respectively, in Fig. 3(a) show that som for sj — 1 and p = 1 is decreased at m > 1 and again 
show the saturation with a value of about 0.87 (0.71) in direct simulation (DMA calculation) 
at m = 20. For < s/ < 1, som show a similar, gradual change as m is increased. It is 
noted that the agreement between the results of DMA calculations and direct simulation is 
good at small m but become worse at larger m. This is due to the adopted NLA in which 
the correlation beyond the nearest layers is neglected. 

It is noted that an increase in the synchrony as m is increased, which is realized for sj = 
and 0.2 in Fig. 3(a), is due to common inputs arising from all-to-all interlayer couplings for 
p = 1 in Eq. (4). In fact, if we set p = for which inputs come only through local couplings 
in Eq. (4), the synchrony is gradually decreased as spikes propagate by effects of random 
noises for all values of s/, as shown in Fig. 3(c). In the intermediate p value, for example, 
for p — 0.4, the synchrony is decreased (increased) compared with that for p = 1 {p — 0), 
as shown in Fig. 3(b). 

Figure 4(a), 4(b) and 4(c) express the m dependence of aom, RMS value of jitters in 
firing times, for p = 1.0, 0.4 and 0, respectively, with various sj values. Although som is 
variable depending on sj and p as shown in Figs. 3(a)-3(c), magnitudes of aom a-re nearly 
independent of m, which shows that spikes propagate with nearly the same dispersion. In 
particular, for p = 0, the m dependence of o"om is almost the same for all sj values, as 
shown in Fig. 4(c). Because of the adopted NLA, the agreement between the results of 
DMA calculations and direct simulation become worse at larger m although both results are 
similar in the qualitative sense. The neuron activity aom defined by Eq. (55) is 0.50 - 0.51 
at m > 1 for all the cases investigated (not shown). 

So far we have adopted values of = 10 and M = 20. It is desirable to perform 
numerical calculations with larger values of N and M for a better understanding of multi- 
layer networks in living brains. Because of a limitation of our computer facility, we have 
performed only DMA calculations for larger value of N and M. Figure 5(a), 5(b) and 5(c) 
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show the m-dependence of som for P — 1-0, 0.4 and 0.0, respectively, with N = 100 and 
M = 40 for various sj values: parameters of f3, Wi, W2 are same as in Fig. 2. When 
comparing Figs. 5(a)-5(c) with 2(a)-2(c), we note similar m dependence in them: the N 
dependence of som will be shortly discussed in Sec. HID. 

B. Effects of p 

As was pointed out in Figs. 3(a)-3(c), the factor of p plays an important role for syn- 
chrony in spike propagation. In order to systematically study the effect of p on the input- 
output relation of M = 20 multilayer, we have calculated So20, Som at m = M = 20, as 
a function of p for various sj values with cx/ = 1, /? = 0.01 and N = 100, whose result is 
shown in Fig. 6. For p — 0, So2o is very small for all sj. When p is increased from 0, So2o 
is linearly increased and it shows an almost saturation at p > 0.5 — 0.6. 

Figure 7(a) depicts the calculated result showing So 20 against sj, the input-output re- 
lation of the correlation for a fixed value of f3 = 0.01 in the multilayer of M = 20 and 

= 100. It is shown that for independent local couplings only (p = 0), <so2o becomes too 
small compared to sj. In contrast, for common all-to-all feedforward couplings only (p = 1), 

50 20 becomes larger than the input correlation for sj < sjc where s/c = 0.54 is the critical 
value below which Sq 20 > Sj. For p = 0.2 and 0.4, the critical value becomes sjc = 0.09 and 
0.33, respectively, which are nearly the same as the experimentally observed value of 0.1-0.3 
[32]- [40]. 

C. Effects of P 

As was shown in Fig. 1(b), noises are detrimental for the synchrony of spikes. This fact 
is realized when we compare som for P = 0.02 shown in Fig. 5(d) with that for /3 = 0.01 in 
Fig. 5(a). The value of som for = 1 at m = 40 is about 0.20 in Fig. 5(b) which is much 
smaller than 0.46 in Fig. 5(a). 

Figure 7(b) expresses sj versus sq 20 when the noise intensity is changed with a fixed 
value of p = 1 for (Tj = 1 and A^ = 100. In the case of P — 0.01, so 20 is larger than sj for 

51 < Sic — 0.54. On the contrary, in the case of (3 — 0.02, the critical value is s/c = 0.18. 
Furthermore, in the case of /? = 0.03, we get sjc — 0.11. Thus sjc is much reduced with 
increasing /3. 

D. Effects of N 

As mentioned above, sq 20 becomes smaller than sj for sj > 0.18 for /3 = 0.02 and 
N — 100. This situation is changed if the size of N is reduced. Figure 8(a) shows som for 
various A^ with sj = 0.4 and f3 = 0.02. In the case of A^ = 10, for example, S020 is 0.50 
which is larger than 0.19 for A" = 100. For A^ = 20, so 20 is nearly the same as s/ = 0.4. 
Figure 8(a) clearly shows that sq 20 is gradually decreased as A^ is increased. 
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E. Effects of wi 



So far we have assumed vanishing intralayer couphngs, Wi, which are now introduced. 
When intralayer couphng wi are positive (excitatory), som is expected to be increased. This 
is confirmed in our calculations shown in Fig. 8(b) depicting som for 1^1=0, 0.05 and 0.1 
with (7/ = 1, Sf = 0.4, f3 = 0.02, W2 = 0.1 and N = 100. On the contrary, if wi is negative 
(inhibitory), it is considered to prevent the propagation of a spike. Actually, Fig. 8(b) shows 
that for w — —0.05, a propagation of a spike is terminated at m = 7 below which som is 
smaller than that for positive Wi. 

F. Effects of 

The interlayer coupling W2 is expected to play also important roles in spike propagation. 
Figure 8(c) shows calculated results when the interlayer coupling W2 is increased from 0.1 
to 0.2, which yields an increase in Som for o'l = 1, Sj = 0.4, Wi = 0, j3 == 0.02 and 
N = 100. On the contrary, our calculation in Fig. 8(c) shows that the negative couplings 
with W2 — —0.1 and -0.2 are not favorable for the spike propagation which is terminated at 
m — 9. 

IV. CONCLUSIONS AND DISCUSSIONS 

We have discussed the spatial correlation while spikes propagate through feedforward 
multilayer. Figures 3(a) and 3(b) suggest that as m is increased, som and aom approach 
fixed values. In order to show this more clearly, we depict, in Fig. 9(a), the (Tom-Som plot in 
which points of (crom, som) are sequentially connected from m = to m = 40 with (5 — 0.01 
and N = 100: note that {<7om, som) for m = stand for input values. For example, in the 
case of o"/ = 1 and s/ = 1.0, the point starts from (o"om, ■5om) = (l-0, 1.0) at m = and ends 
with (0.58, 0.45) at m = 40. In contrast, in the case of cx/ = 1 and sj = 0.0, the point starts 
from (1.0, 0.0) at m = and ends with (0.49, 0.22) at m = 40. In the case oi aj — and 
Si = 0.0, the point varies from (0.0, 0.0) at m = to (0.48, 0.21) at m = 40. These show 
that a fixed point may be about {crooo, sooo) ~ (0.54,0.38), as shown by the cross in Fig. 
9(a). 

Figure 9(b) show a similar aom-som piot for a larger P = 0.02. In the case of a/ = 1 
and Si — 1.0, for example, the point starts from (1.0, 1.0) at m = and ends with (0.95, 
0.22) at m = 40. In the case oi ai — 1 and s/ = 0.0, the point changes from (10.0, 1.0) at 
m = to (0.92, 0.16) at m = 40. Results for ai = and s/ = and for ai — 1 (and 2) with 
< s/ < 1 show that (ciooo, sooo) ~ (0.93, 0.18) for p = 0.02. 

Including calculated results of the neuron activity aom [Eq. (55)], which becomes 
O'Om — 0.50 — 0.51 at m > 1 for all the cases investigated (not shown), we may say that 
all curves starting from different initial values of o"/ and s/ converge to the fixed point of 
{aooo, cooo, sood) in the three-dimensional space spanned by aom, CFOm and som- The fixed 
point is determined by the parameters characterizing the multilayer architecture such as 
/5, p, wi, W2 and A^, but independently of the parameters of o"/ and si for input signals. 
Our conclusion supplements the result of Diesmann et al. [18] who have shown that in the 
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propagating regime, the number of firing neurons and RMS of firing times in a pulse packet 
converge to fixed-point values. Our calculation has shown that while a pulse packet propa- 
gates with an almost constant dispersion (RMS), the spatial correlation within the packet 
may change, and in a deep layer, it saturates at the value determined by the parameters 
depending on the multilayer. 

The dependence of calculated fixed-point values of sooo and aooo on /3, p, N, wi and W2 
are summarized as follows. 

(1) o"ooo is increased as increasing (3, but decreased as increasing p, N, wi or W2- 

(2) sooo is increased as increasing p, wi or W2, but decreased as increasing (3 or N. 

We have tried to elucidate this property by an analysis using DMA. Because the fixed points 
do not depend on aj and sj, we consider the case of aj — sr — 0. In the case of wi — W2 — 0, 
Eqs. (28)- (52) yield 

oc /3^ (61) 

pT « I' (62) 

for m — > oo where they are independent of m. When wi and W2 are small, Eqs. (30), (33), 
(37), (39), (45) and (48) yield following equations given as series of wi and W2: 

7i,i = lim^^oo 7™i™ = C/3^(l - aiWi - 02^2) +dw2\pp[^i + (1 -p) (63) 
= lim^^oo pT,r = (^ ~ ^1^1 ~ ^2«^2) + d W2p'i_i, (64) 

7^1 = lim^^oo 7™f = ew2[p Pi,i + {I - p) 71,1], (65) 

' !• 171—1,171 tCC\ 

= lim^^oo Pi,i =ew2 pi,i, (66) 

where expansion coefficients of ai, 02, &i, &2 c, d and e are obtainable from Eqs. (28)-(52) 
in principle although their explicit forms are not necessary for our qualitative discussion. 
Solving Eqs. (63)- (66) for 71^1 and which are substituted to Eqs. (57) and (60), we get 

crooc oc j (1 - ^(aiwi + 02^2) + ^dewl [(1 - pf + ^p{2 - p)]^ , (67) 

(ai-bi\ (a2-h2\ (dep{2-p)\ ^ 
.000 - [j^ j -1 + [-j^ j W2 + j ^2, (68) 

where only relevant terms are retained. Expressions given by Eqs. (67) and (68) may 
account for all the dependence of aooo and sqoo raised in items (1) and (2), except the N 
dependence of aooo and the f3 dependence of sqoc- The former may be explained if Oi or 02 
is an increasing function of N, and the latter if ai or a2 (61 or 62) is a decreasing (increasing) 
function of (3. 

Although numerical calculations reported in Sec. Ill have been made for single spike 
inputs, we may easily apply our DMA theory to the case of spike-train inputs given by [see 
Eq. (8)] 

I^"\t) =6,niuJ2'^it-tije) (69) 
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where tjjg stands for the £th input time to neuron j. Raster in Fig. 10(a) shows firings of 
neurons for an apphed Poisson spike train with the average interspike interval (ISI) of 100 
{p = 0.02 ai = 1, Si = 0, p = 1, Wi = 0, W2 = 0.1, M = 20 and iV = 10) calculated by 
direct simulation (a single run). In contrast, raster in Fig. 10(b) expresses ^Zo&a^ firings on 
layer m for which the firing time is defined by 

tomg = {t\i^Tit)=0;i:iT>o}. (70) 

where n'^it) denotes the averaged voltage variable on the layer m [Eq. (17)]. Time courses 
of fi^{t) are plotted in Fig. 10(c), in which the uppermost frame shows an applied Poisson 
spike train, I^^\t). The resuh of ^^^(t) of DMA (sohd curves) is in good agreement with that 
of direct simulations (dashed curves). Figures 10 (a)- (c) show that spikes propagate from 
m = 1 to m = 10 of the M = 20 multilayer. Spurious firings due to added noises, which 
is reahzed in Fig. 10(a), vanish in global firings shown in Fig. 10(b) by the averaging over 
neuron ensembles, which expresses the population effect [13] [14]. Comparing fi^{t) with 
an applied spike train of /'^^^(t), for example, at t ~ 300 — 600, we note that when ISI of 
input spikes is shorter than about 55, FN neuron cannot respond because of its refractory 
period, which is realized also in HH neuron [51]. It is noted that although the correlation 
for spike-train inputs develops while spikes propagate through the multilayer, as for single 
spike inputs, means and variances of their ISI remain almost the constant. 

To summarize, we have studied the spatial correlation during spike's propagation through 
multilayers to show 

(a) the input correlation of s/ may propagate through the network, yielding som ~ at the 
end layer of m = M ~ 10 — 20 with the observed magnitude of sqm ~ 0.1 — 0.3 [32]- [40], 
when model parameters are appropriate, and 

(b) in a long multilayer, the correlation of the deep layer converges to a fixed-point value of 
(cooo, sooo) which depends on the parameters characterizing the multilayer architecture but 
is independent of the input correlation. 

The item (1) implies that spikes in multilayers with physiologically reasonable size of m = 
M ~ 10 — 20 may carry information encoded in the spatial correlation of firing times. The 
item (b) is similar to the result obtained in Rcf. [23], where the spike rate in a deep layer of 
a balanced synfire chain is shown to be independent of the input rate. 

Finally we would hke to point out the efficiency of our DMA. Direct simulations with 
100 trials for a multilayer of M = 20 and N — 10 with a set of parameters (Figs. 1 and 
2), required the computation time of 51 minutes by using 1.8 GHz CPU PC, while a DMA 
calculation needs only 6 s, which is about 500 times faster than simulations. 
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APPENDIX A: DERIVATION OF EQS. (41)-(44) 

From Eqs. (l)-(5), we get DEs for the deviations of 5xmj and 5ymj of a neuron j (=1 to 
N) in the layer m {—1 to M), given by (see Appendix A in I) 



dt 
dSymj 
dt 

with 



= frSxmj + mSxl^ - 7m) + fr^^j - cSymj + Sli^f + 5ltf + 5lt] + Uj, (Al) 
= hSxmj - ddymj, {j e m) (A2) 

= (f^) i: (A3) 



Sli^f = (1 - 6ml) W2gT[^J2^Xm-lk + {I - p)SXm-lji (A4) 



= -5^1 u H 6tij, (A5) 

where H = dH(t)/d,t. Wc have taken into account up to third-order terms in dxmj which 
play an important role in stabilizing Des [24], while only a linear-order term is included in 
the coupling term in Eqs. (A3) and (A4). DEs for the variances and covariances are given 
by 

^ = ^ E[2 < Sxm, (^) > 6^,6,,+ < Sy^j (^) + Sx^, (^) > S^,S,, 
j 

+ 2<5ymj{^)> SM (A6) 
^-^.J:U^< S^nj (^) > 5^,5,,+ < 5ym, (^) + Sx^, (^) > 5.^5,, 

+ < Synj (^) + SXmk (^) > 5.25ai + 2 < (^) > (A7) 
dt dt dt 

Substituting Eqs. (A1)-(A5) to Eqs. (A6) and (A7), we get DEs for 7^^'™ and p^''^" {n, A = 
1,2). In the process of these calculations, we get new correlation functions of Pnit) and 
Rnit) defined by 

Pnit) = E(< ^^ijit) Stij > (5,1+ < Syijit) Stij > S,2), (A8) 
j 

^«(^) ^ -^1212i< ^^uW ^^ik > < ^yijit) stik > 5«2), (A9) 

whose equations of motions are given by Eqs. (41)-(44). 
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FIGURES 



FIG. 1. Raster showing firings of neurons on the first ten layers in a multilayer oi N = 10 and 

M = 20 for (a) /3 = 0.01 and (b) (3 = 0.02 with aj = 1, sj = 0, p = 1, wi = 0, W2 = 0.1, M = 20 
and N = 10, calculated by a direct simulation (a single trial). The vertical scale expresses the 
neuron index k defined hy k = 10(m — 1) + j where 1 < m < 10 and 1 < j < 10. 

FIG. 2. Time courses of (a) /x™ and (b) Som in 1 < m < 10 of a multilayer with iV = 10 and 
M = 20 for ai = 1, SI = 0, p = I, P = 0.01, wi = and W2 = 0.1, calculated by DMA (solid 
curves) and direct simulations of 100 trials (dashed curves). 

FIG. 3. som for (a) p = 1.0, (b) p = 0.4 and (c) p = 0.0, with a/ = 1, /? = 0.01, wi = 0, 
W2 = 0.1, A'^ = 10 and various s/: sj = 1.0 (squares), 0.8 (circles), 0.6 (diamonds), 0.4 (inverted 
triangles), 0.2 (triangle) and (large circles), calculated by direct simulations (open marks) and 
DMA (filled marks): 

FIG. 4. aom for (a) p = 1.0, (b) p = 0.4 and (c) p = 0.0, with ai = 1, p = 0.01, wi = 0, 
W2 = 0.1, N = 10 and various sj: sj = 1.0 (squares), 0.8 (circles), 0.6 (diamonds), 0.4 (inverted 
triangles), 0.2 (triangle) and (large circles), calculated by direct simulations (open marks) and 
DMA (filled marks): 

FIG. 5. som for (a) /? = 0.01 and p = 1.0, (b) f3 = 0.01 and p = 0.4, (c) /3 = 0.01 and p = 0.0, 
and (d) /? = 0.02 and p = 1, with aj = 1, wi = 0, W2 = 0.1, N = 100 and various s/ calculated by 
DMA: sj = 1.0 (squares), 0.8 (circles), 0.6 (diamonds), 0.4 (inverted triangles), 0.2 (triangle) and 
(large circles). 

FIG. 6. The p dependence of so 20 for various sj with aj = 1, (3 = 0.01, wi = 0, W2 = 0.1 
and N = 100 calculated by DMA: sj = 0.80 (circles), 0.6 (diamonds), 0.4 (inverted triangles), 0.2 
(triangle) and (large circles). 

FIG. 7. (a) 5(9 20 against sj for /3 = 0.01 with various p: p = 1 (circles), 0.6 (squares), 0.4 
(inverted triangles), 0.2 (triangles) and (diamonds), (b) so 20 against s/ for p = 1 with various 
(5: j3 = 0.01 (circles), 0.02 (triangles) and 0.03 (squares), (a) and (b) are calculated by DMA with 
crj = l,wi= 0, W2 = 0.1 and N = 100. 

FIG. 8. (a) som for different N, (b) for different wi, and (c) for different W2, with aj = 1, 
si = 0.4, p = 0.4, p = 0.02, calculated by DMA. 

FIG. 9. aom against som for (a) f3 = 0.01 and (b) /3 = 0.02, with aj = 1, p = 1 wi = 0, 
W2 = 0.1, M = 40 and N = 100 calculated by DMA. Arrows denote the direction of increasing m. 
All points starting from {aoo, soo) converge to the fixed-point marked by the cross (see text). 
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FIG. 10. (a) Raster of firings of individual neurons, (b) raster of global firings averaged on 
each layer in a multilayer calculated by a direct simulation (a single trial), and (c) time courses 
of calculated by DMA (solid curves) and simulations (dashed curves) for applied Poisson spike 
inputs /^^^ shown at the uppermost frame of (c): /3 = 0.02 aj = 1, sr = 0, p = 1, wi = 0, 
1^2 = 0.1, M = 20 and N = 10. The vertical scale of (a) expresses the neuron index k defined by 
k = 10(m — 1) + j where 1 < m < 10 and 1 < j < 10 (see text). 
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